library(MASS)
library(Zelig)
library(MNP)
#library(scatterplot3d)
#rm(list=ls())

load("home\\rchaudoin\\chaudoin\\Interact_031412_long.RData")

#############################################################################
####  Predictions: Unemployment, Election Years
#############################################################################

### This constructs a marix se1a, which is the "newdata" for predicted probabilities-- PE = 1
se1a<-WTOdata[0:121,]
se1a$PE1_UE<-as.matrix(seq(4.5,5.7,by=.01))
se1a$Unemp_6mo<-as.matrix(seq(4.5,5.7,by=.01))
se1a$PE1<-as.matrix(rep(1,121), nrow=121, ncol=1)
se1a$Exports_6mo<-as.matrix(rep(mean(WTOdata$Exports_6mo),121), nrow=121, ncol=1)
se1a$Imports_6mo<-as.matrix(rep(mean(WTOdata$Imports_6mo),121), nrow=121, ncol=1)
se1a$Plaintiff_PCGDP<-as.matrix(rep(mean(WTOdata$Plaintiff_PCGDP,na.rm=TRUE),121), nrow=121, ncol=1)
se1a$Plaintiff_UE<-as.matrix(rep(mean(WTOdata$Plaintiff_UE,na.rm=TRUE),121), nrow=121, ncol=1)
se1a$PlEP1<-as.matrix(rep(1,121), nrow=121, ncol=1)
se1a$PlEP1_PlUE<-as.matrix(rep(mean(WTOdata$Plaintiff_UE,na.rm=TRUE),121), nrow=121, ncol=1)
se1a$Age<-as.matrix(rep(mean(WTOdata$Age),121), nrow=121, ncol=1)
se1a$Age2<-as.matrix(rep(mean(WTOdata$Age2),121), nrow=121, ncol=1)
se1a$Month<-as.matrix(rep(mean(WTOdata$Month),121), nrow=121, ncol=1)
se1a$Month2<-as.matrix(rep(mean(WTOdata$Month2),121), nrow=121, ncol=1)


# This constructs se1b, which is the "newdata" for predicted probabilities-- PE = 0
se1b<-se1a
se1b$PE1_UE<-as.matrix(rep(0,121),nrow=121,ncol=1)
se1b$PE1<-as.matrix(rep(0,121),nrow=121,ncol=1)


###########
# Predicted Probabilities Loops: UE/PE1-0
###########

### First, PE = 1

# Runs the predict command once... I'll cbind the predictions together
pred_m1 <-predict(m3, newdata = se1a, type = "prob")
pred_m1_0 <- pred_m1$p[,1]
pred_m1_3 <- pred_m1$p[,2]
pred_m1_6 <- pred_m1$p[,3]

# Loops over the predict command z times, creates matrix with z columns, each column a new prediction
q<-1
	while (q<11) {
	z<-1
	while (z<99) {
		print (z)
		pred_m1_junk <-predict(m3, newdata = se1a, type = "prob")
		pred_m1_0 <- cbind(pred_m1_0, pred_m1_junk$p[,1])
		pred_m1_3 <- cbind(pred_m1_3, pred_m1_junk$p[,2])
		pred_m1_6 <- cbind(pred_m1_6, pred_m1_junk$p[,3])
		z <- z+1
		}
	q <- q+1
	}

# Making a 121x4 matrix that will have the values of UE in c1 and means of the pred probs in c2-c4
se_m1_ue<-as.matrix(seq(4.5,5.7,by=.01))

# For now, the columns for predicted probabilities are just the sequence to keep tabs/dimensions straight
pred_m1_mean0_pe1<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_mean3_pe1<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_mean6_pe1<-as.matrix(seq(4.5,5.7,by=.01))
# Upper 95%
pred_m1_ub0_pe1<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_ub3_pe1<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_ub6_pe1<-as.matrix(seq(4.5,5.7,by=.01))
# Lower 95%
pred_m1_lb0_pe1<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_lb3_pe1<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_lb6_pe1<-as.matrix(seq(4.5,5.7,by=.01))


# Takes the mean of each row of predictions, UB and LB
z<-1
while (z<122) {
	# Means
	pred_m1_mean0_pe1[z,] <- mean(pred_m1_0[z,])
	pred_m1_mean3_pe1[z,] <- mean(pred_m1_3[z,])
	pred_m1_mean6_pe1[z,] <- mean(pred_m1_6[z,])
	# UBs
	pred_m1_ub0_pe1[z,] <- quantile(pred_m1_0[z,], probs=.95)
	pred_m1_ub3_pe1[z,] <- quantile(pred_m1_3[z,], probs=.95)
	pred_m1_ub6_pe1[z,] <- quantile(pred_m1_6[z,], probs=.95)
	# LBs
	pred_m1_lb0_pe1[z,] <- quantile(pred_m1_0[z,], probs=.05)
	pred_m1_lb3_pe1[z,] <- quantile(pred_m1_3[z,], probs=.05)
	pred_m1_lb6_pe1[z,] <- quantile(pred_m1_6[z,], probs=.05)
	z<-z+1
	}

###
# PE = 1, prediction data
# Creates a matrix with UE and the means of the predictions, UBs, and LBs
se_m1_full_pe1<-cbind(se_m1_ue, pred_m1_mean0_pe1, pred_m1_mean3_pe1, pred_m1_mean6_pe1, pred_m1_ub0_pe1, pred_m1_ub3_pe1, pred_m1_ub6_pe1, pred_m1_lb0_pe1, pred_m1_lb3_pe1, pred_m1_lb6_pe1)

# Test plots
#plot(se_m1_full_pe1[,1],se_m1_full_pe1[,4], type="l", ylim=c(0,0.005))
#lines(se_m1_full_pe1[,1],se_m1_full_pe1[,7], type="l", col=3)
#lines(se_m1_full_pe1[,1],se_m1_full_pe1[,10], type="l", col=4)


##################################
# Now, PE = 0
####
# Runs the predict command once... I'll cbind the predictions together
pred_m1_pe0 <-predict(m3, newdata = se1b, type = "prob")
pred_m1_0_pe0 <- pred_m1_pe0$p[,1]
pred_m1_3_pe0 <- pred_m1_pe0$p[,2]
pred_m1_6_pe0 <- pred_m1_pe0$p[,3]

# Loops over the predict command z times, creates matrix with z columns, each column a new prediction
q<-1
	while (q<11) {
		z<-1
		while (z<99) {
			print (z)
			pred_m1_junk_pe0 <-predict(m3, newdata = se1b, type = "prob")
			pred_m1_0_pe0 <- cbind(pred_m1_0_pe0, pred_m1_junk_pe0$p[,1])
			pred_m1_3_pe0 <- cbind(pred_m1_3_pe0, pred_m1_junk_pe0$p[,2])
			pred_m1_6_pe0 <- cbind(pred_m1_6_pe0, pred_m1_junk_pe0$p[,3])
			z <- z+1
		}
	q <- q+1
	}

# Making a 121x4 matrix that will have the values of UE in c1 and means of the pred probs in c2-c4
se_m1_ue<-as.matrix(seq(4.5,5.7,by=.01))

# For now, the columns for predicted probabilities are just the sequence to keep tabs/dimensions straight
pred_m1_mean0_pe0<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_mean3_pe0<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_mean6_pe0<-as.matrix(seq(4.5,5.7,by=.01))
# Upper 95%
pred_m1_ub0_pe0<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_ub3_pe0<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_ub6_pe0<-as.matrix(seq(4.5,5.7,by=.01))
# Lower 95%
pred_m1_lb0_pe0<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_lb3_pe0<-as.matrix(seq(4.5,5.7,by=.01))
pred_m1_lb6_pe0<-as.matrix(seq(4.5,5.7,by=.01))


# Takes the mean of each row of predictions, UB and LB
z<-1
while (z<122) {
	# Means
	pred_m1_mean0_pe0[z,] <- mean(pred_m1_0_pe0[z,])
	pred_m1_mean3_pe0[z,] <- mean(pred_m1_3_pe0[z,])
	pred_m1_mean6_pe0[z,] <- mean(pred_m1_6_pe0[z,])
	# UBs
	pred_m1_ub0_pe0[z,] <- quantile(pred_m1_0_pe0[z,], probs=.95)
	pred_m1_ub3_pe0[z,] <- quantile(pred_m1_3_pe0[z,], probs=.95)
	pred_m1_ub6_pe0[z,] <- quantile(pred_m1_6_pe0[z,], probs=.95)
	# LBs
	pred_m1_lb0_pe0[z,] <- quantile(pred_m1_0_pe0[z,], probs=.05)
	pred_m1_lb3_pe0[z,] <- quantile(pred_m1_3_pe0[z,], probs=.05)
	pred_m1_lb6_pe0[z,] <- quantile(pred_m1_6_pe0[z,], probs=.05)
	z<-z+1
	}

###
# PE = 0, prediction data
# Creates a matrix with UE and the means of the predictions, UBs, and LBs
se_m1_full_pe0<-cbind(se_m1_ue, pred_m1_mean0_pe0, pred_m1_mean3_pe0, pred_m1_mean6_pe0, pred_m1_ub0_pe0, pred_m1_ub3_pe0, pred_m1_ub6_pe0, pred_m1_lb0_pe0, pred_m1_lb3_pe0, pred_m1_lb6_pe0)


save.image("home\\rchaudoin\\chaudoin\\Interact_2013_03_05_long_PEUE")











##########################################################################
# Effects of Exports
##########################################################################


############################
# Making newdata for EX
############################

se1c<-WTOdata[0:121,]
se1c$Exports_6mo<-as.matrix(seq(1,4.96,by=.033))
se1c$PE1_UE<-as.matrix(rep(mean(WTOdata$Unemp_6mo),121), nrow=121, ncol=1)
se1c$Unemp_6mo<-as.matrix(rep(mean(WTOdata$Unemp_6mo),121), nrow=121, ncol=1)
se1c$PE1<-as.matrix(rep(1,121), nrow=121, ncol=1)
se1c$Imports_6mo<-as.matrix(rep(mean(WTOdata$Imports_6mo),121), nrow=121, ncol=1)
se1c$Plaintiff_PCGDP<-as.matrix(rep(mean(WTOdata$Plaintiff_PCGDP,na.rm=TRUE),121), nrow=121, ncol=1)
se1c$Plaintiff_UE<-as.matrix(rep(mean(WTOdata$Plaintiff_UE,na.rm=TRUE),121), nrow=121, ncol=1)
se1c$PlEP1<-as.matrix(rep(1,121), nrow=121, ncol=1)
se1c$PlEP1_PlUE<-as.matrix(rep(mean(WTOdata$Plaintiff_UE,na.rm=TRUE),121), nrow=121, ncol=1)
se1c$Age<-as.matrix(rep(mean(WTOdata$Age),121), nrow=121, ncol=1)
se1c$Age2<-as.matrix(rep(mean(WTOdata$Age2),121), nrow=121, ncol=1)
se1c$Month<-as.matrix(rep(mean(WTOdata$Month),121), nrow=121, ncol=1)
se1c$Month2<-as.matrix(rep(mean(WTOdata$Month2),121), nrow=121, ncol=1)


###
# Predicted Probabilities Loops: Exports
###

# Runs the predict command once... I'll cbind the predictions together
pred_m3_ex <-predict(m3, newdata = se1c, type = "prob")
pred_m3_ex0 <- pred_m3_ex$p[,1]
pred_m3_ex3 <- pred_m3_ex$p[,2]
pred_m3_ex6 <- pred_m3_ex$p[,3]

# Loops over the predict command z times, creates matrix with z columns, each column a new prediction
z<-1
while (z<300) {
	print (z)
	pred_m3_junk <-predict(m3, newdata = se1c, type = "prob")
	pred_m3_ex0 <- cbind(pred_m3_ex0, pred_m3_junk$p[,1])
	pred_m3_ex3 <- cbind(pred_m3_ex3, pred_m3_junk$p[,2])
	pred_m3_ex6 <- cbind(pred_m3_ex6, pred_m3_junk$p[,3])
	z <- z+1
	}

# Making a 121x4 matrix that will have the values of EX in c1 and means of the pred probs in c2-c4
se1c_ex<-as.matrix(seq(1,4.96,by=.033))

# For now, the columns for predicted probabilities are just the sequence to keep tabs/dimensions straight
pred_m3_mean0_ex<-as.matrix(seq(1,4.96,by=.033))
pred_m3_mean3_ex<-as.matrix(seq(1,4.96,by=.033))
pred_m3_mean6_ex<-as.matrix(seq(1,4.96,by=.033))
# Upper 95%
pred_m3_ub0_ex<-as.matrix(seq(1,4.96,by=.033))
pred_m3_ub3_ex<-as.matrix(seq(1,4.96,by=.033))
pred_m3_ub6_ex<-as.matrix(seq(1,4.96,by=.033))
# Lower 95%
pred_m3_lb0_ex<-as.matrix(seq(1,4.96,by=.033))
pred_m3_lb3_ex<-as.matrix(seq(1,4.96,by=.033))
pred_m3_lb6_ex<-as.matrix(seq(1,4.96,by=.033))


# Takes the mean of each row of predictions- mean prediction for each outcome
z<-1
while (z<122) {
	# Means
	pred_m3_mean0_ex[z,] <- mean(pred_m3_ex0[z,])
	pred_m3_mean3_ex[z,] <- mean(pred_m3_ex3[z,])
	pred_m3_mean6_ex[z,] <- mean(pred_m3_ex6[z,])
	# UBs
	pred_m3_ub0_ex[z,] <- quantile(pred_m3_ex0[z,], probs=.95)
	pred_m3_ub3_ex[z,] <- quantile(pred_m3_ex3[z,], probs=.95)
	pred_m3_ub6_ex[z,] <- quantile(pred_m3_ex6[z,], probs=.95)
	# LBs
	pred_m3_lb0_ex[z,] <- quantile(pred_m3_ex0[z,], probs=.05)
	pred_m3_lb3_ex[z,] <- quantile(pred_m3_ex3[z,], probs=.05)
	pred_m3_lb6_ex[z,] <- quantile(pred_m3_ex6[z,], probs=.05)
	z<-z+1
	}

# Creates a matrix with EX and the means of the predictions, ub, lb
se_m3_full_ex<-cbind(se1c_ex, pred_m3_mean0_ex, pred_m3_mean3_ex, pred_m3_mean6_ex, pred_m3_ub0_ex, pred_m3_ub3_ex, pred_m3_ub6_ex, pred_m3_lb0_ex, pred_m3_lb3_ex, pred_m3_lb6_ex)

# Saving
save.image("home\\rchaudoin\\chaudoin\\Interact_2013_03_05_long_EX")










##########################################################################
# Effects of Imports
##########################################################################

############################
# Making newdata for IM
############################

# Effects of Imports
se1d<-WTOdata[0:121,]
se1d$Imports_6mo<-as.matrix(seq(0.9,7.62,by=.056))
se1d$PE1_UE<-as.matrix(rep(mean(WTOdata$Unemp_6mo),121), nrow=121, ncol=1)
se1d$Unemp_6mo<-as.matrix(rep(mean(WTOdata$Unemp_6mo),121), nrow=121, ncol=1)
se1d$PE1<-as.matrix(rep(1,121), nrow=121, ncol=1)
se1d$Exports_6mo<-as.matrix(rep(mean(WTOdata$Exports_6mo),121), nrow=121, ncol=1)
se1d$Plaintiff_PCGDP<-as.matrix(rep(mean(WTOdata$Plaintiff_PCGDP,na.rm=TRUE),121), nrow=121, ncol=1)
se1d$Plaintiff_UE<-as.matrix(rep(mean(WTOdata$Plaintiff_UE,na.rm=TRUE),121), nrow=121, ncol=1)
se1d$PlEP1<-as.matrix(rep(1,121), nrow=121, ncol=1)
se1d$PlEP1_PlUE<-as.matrix(rep(mean(WTOdata$Plaintiff_UE,na.rm=TRUE),121), nrow=121, ncol=1)
se1d$Age<-as.matrix(rep(mean(WTOdata$Age),121), nrow=121, ncol=1)
se1d$Age2<-as.matrix(rep(mean(WTOdata$Age2),121), nrow=121, ncol=1)
se1d$Month<-as.matrix(rep(mean(WTOdata$Month),121), nrow=121, ncol=1)
se1d$Month2<-as.matrix(rep(mean(WTOdata$Month2),121), nrow=121, ncol=1)

###
# Predicted Probabilities Loops: Imports
###

# Runs the predict command once... I'll cbind the predictions together
pred_m3_im <-predict(m3, newdata = se1d, type = "prob")
pred_m3_im0 <- pred_m3_im$p[,1]
pred_m3_im3 <- pred_m3_im$p[,2]
pred_m3_im6 <- pred_m3_im$p[,3]

# Loops over the predict command z times, creates matrix with z columns, each column a new prediction
z<-1
while (z<300) {
	print (z)
	pred_m3_junk <-predict(m3, newdata = se1d, type = "prob")
	pred_m3_im0 <- cbind(pred_m3_im0, pred_m3_junk$p[,1])
	pred_m3_im3 <- cbind(pred_m3_im3, pred_m3_junk$p[,2])
	pred_m3_im6 <- cbind(pred_m3_im6, pred_m3_junk$p[,3])
	z <- z+1
	}

# Making a 121x4 matrix that will have the values of EX in c1 and means of the pred probs in c2-c4
se1d_im<-as.matrix(seq(0.9,7.62,by=.056))

# For now, the columns for predicted probabilities are just the sequence to keep tabs/dimensions straight
pred_m3_mean0_im<-as.matrix(seq(0.9,7.62,by=.056))
pred_m3_mean3_im<-as.matrix(seq(0.9,7.62,by=.056))
pred_m3_mean6_im<-as.matrix(seq(0.9,7.62,by=.056))
# Upper 95%
pred_m3_ub0_im<-as.matrix(seq(0.9,7.62,by=.056))
pred_m3_ub3_im<-as.matrix(seq(0.9,7.62,by=.056))
pred_m3_ub6_im<-as.matrix(seq(0.9,7.62,by=.056))
# Lower 95%
pred_m3_lb0_im<-as.matrix(seq(0.9,7.62,by=.056))
pred_m3_lb3_im<-as.matrix(seq(0.9,7.62,by=.056))
pred_m3_lb6_im<-as.matrix(seq(0.9,7.62,by=.056))


# Takes the mean of each row of predictions- mean prediction for each outcome
z<-1
while (z<122) {
	# Means
	pred_m3_mean0_im[z,] <- mean(pred_m3_im0[z,])
	pred_m3_mean3_im[z,] <- mean(pred_m3_im3[z,])
	pred_m3_mean6_im[z,] <- mean(pred_m3_im6[z,])
	# UBs
	pred_m3_ub0_im[z,] <- quantile(pred_m3_im0[z,], probs=.95)
	pred_m3_ub3_im[z,] <- quantile(pred_m3_im3[z,], probs=.95)
	pred_m3_ub6_im[z,] <- quantile(pred_m3_im6[z,], probs=.95)
	# LBs
	pred_m3_lb0_im[z,] <- quantile(pred_m3_im0[z,], probs=.05)
	pred_m3_lb3_im[z,] <- quantile(pred_m3_im3[z,], probs=.05)
	pred_m3_lb6_im[z,] <- quantile(pred_m3_im6[z,], probs=.05)
	z<-z+1
	}

# Creates a matrix with EX and the means of the predictions, ub, lb
se_m3_full_im<-cbind(se1d_im, pred_m3_mean0_im, pred_m3_mean3_im, pred_m3_mean6_im, pred_m3_ub0_im, pred_m3_ub3_im, pred_m3_ub6_im, pred_m3_lb0_im, pred_m3_lb3_im, pred_m3_lb6_im)

# Saving
save.image("home\\rchaudoin\\chaudoin\\Interact_2013_03_05_long_IM")